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ABSTRACT 


Noise  generated  by  underwater  sources  is  typically  characterized  in  terms  of  a  far-field,  plane- 
wave  equivalent  source  level  based  on  measurements  made  in  a  free- field  environment  such  as  a 
deep  water  test  range.  Measurements  made  in  a  harbor  environment,  where  multiple  reflections, 
high  background  noise  and  short  propagation  paths  are  typical,  violates  these  conditions.  With 
careful  analysis  it  may  be  possible  to  arrive  at  a  valid  estimate  of  source  location  and  source 
level  from  such  measurements. 

This  paper  reports  results  from  a  test  conducted  at  the  US  Navy  Acoustic  Research  Detachment 
in  Bayview,  Idaho  during  the  summers  of  2010  and  2011.  A  line  array  of  omnidirectional 
hydrophones  was  deployed  from  a  barge.  A  series  of  test  signals  was  transmitted  in  both  deep 
and  shallow  water  using  calibrated  acoustic  sources  to  evaluate  the  effectiveness  of  post¬ 
processing  techniques,  as  well  as  line  array  beamforming,  in  minimizing  reflected  path 
contributions  and  improving  signal-to-noise  ratio.  A  method  of  estimating  the  location  of  the 
sources  while  taking  into  account  a  real,  non  linear  array  based  on  these  measurements  is 
presented. 
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LONG-TERM  GOALS  &  OBJECTIVES 


Source  characterization  measurements  are  best  done  in  a  deep  water  environment  where  far- field 
and  free- field  assumptions  can  be  made.  In  this  case,  acoustic  waves  arriving  at  the  array  are 
approximately  planar,  and  plane  wave  beamforming  can  be  utilized.  This  is  the  only  method 
covered  in  most  array  signal  processing  books.  Plane  wave  beamforming  requires  that  the 
source  is  a  radiating  simple  source  and  the  wavefronts  received  by  the  array  are  planar.  The  far- 
field  is  often  defined  as  starting  at  [1,  2] 


2  ]} 
r  =  — 

A  (i) 

where  r  is  the  distance  from  the  center  of  the  array  to  the  source,  L  is  the  largest  dimension  of  the 
array,  and  X  is  the  acoustic  wavelength.  This  is  the  definition  that  is  used  throughout  this  paper. 

When  the  source  is  closer  to  the  array,  wavefront  curvature  must  be  taken  into  consideration. 
Different  methods  of  accounting  for  this  complexity  have  been  established  [1,  2,  3,  4,  5,  6], 
however  few  have  taken  into  account  the  possibility  of  a  curved  one-dimensional  array  and 
possible  error  in  the  location  of  the  array  [1,2,  3],  and  none  have  used  the  near  field 
compensation  method  of  beamforming  to  account  for  this.  In  this  paper,  a  means  of  source 
localization  is  developed  which  would  work  in  a  shallow,  multipath  environment  with  the  source 
in  the  near- field  of  the  array. 

APPROACH 

This  work  will  present  data  acquired  at  the  Navy’s  Acoustics  Research  Detachment  at  Lake  Pend 
Oreille  in  Bayview,  Idaho.  A  14  element  array  was  deployed  in  a  45  foot  deep  water  column  as 
well  as  in  deep  water  (approximately  1200  feet).  Two  sources  were  also  deployed  and  a  variety 
of  signals  were  transmitted  from  each  source.  The  data  is  processed  using  the  near  field 
compensation  method  of  beamforming  assuming  a  vertical  array  and  compared  to  results  found 
after  using  a  cross-correlation  to  determine  the  actual  shape  of  the  array. 

The  deep  water  test  site  is  examined  first  with  the  array  in  the  far  field  of  the  source.  This  is  the 
simplest  situation,  and  evaluation  using  delay  and  sum  beamforming  in  the  frequency  domain  is 
done  to  ensure  that  the  algorithm  is  working  properly.  This  plane-wave  beamformer  is  proven  to 
work  well  in  this  free  field  environment  when  the  array  is  in  the  far  field  of  the  source.  The 
near  field  compensation  method  of  beamforming  is  then  used  for  this  case.  It  is  shown  that  there  is 
little  to  no  difference  in  the  results  for  both  methods.  Showing  the  near  field  beamformer 
matches  the  plane  wave  beamformer  in  the  far  field  case  proves  that  the  near  field  beamformer  is 
working  properly. 


Following  these  analyses,  the  deep  water  test  site  is  once  again  utilized,  but  with  the  array  in  the 
near  field  of  the  source.  This  allows  for  examination  of  the  effect  of  the  array  being  in  the  near 
field  without  the  added  complication  of  multiple  reflections.  First  the  plane-wave  beamformer 
is  used,  and  is  shown  to  give  inaccurate  results.  The  near  field  algorithm  is  then  used,  giving 
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accurate  results  and  proving  that  it  works  in  near  field  conditions.  There  is  a  small  error  in  the 
location  estimate,  so  a  measurement  of  the  array  element  locations  is  made  to  increase  accuracy. 

The  algorithm  is  then  adjusted  for  the  different  element  locations  and  a  new  set  of  results  is 
found  that  increases  accuracy. 

Finally,  the  array  and  sources  are  moved  into  a  shallow  water,  harbor  environment.  This  is  a 
complex  environment  with  many  reflections  and  the  array  is  in  the  near  field  of  the  source. 
The  original  near  field  algorithm  is  tested  as  well  as  the  algorithm  that  has  been  adjusted  for 
curvature  of  the  array.  In  the  shallow  water  environment  in  particular  the  curvature  is  very 
noticeable  and  making  the  adjustments  for  it  shows  to  have  increased  accuracy. 

OBJECTIVES  (N/A) 

WORK  COMPLETED 

Multiple  methods  of  beam  forming,  both  far  field  and  near  field,  are  discussed  as  a  background  as 
well  as  to  present  methods  used  in  this  paper.  In  the  far  field  regime,  standard  time  delay 
beamforming  is  explored  in  both  the  time  and  frequency  domains.  Two  other  forms  of  far  field 
beamforming  that  are  presented  in  [3]  are  also  discussed,  Bartlett  beamforming  and  an  adaptive 
method  called  the  Minimum  Variance  distortionless  processor  (MV). 

Following  these  far  field  approaches,  several  near  field  algorithms  are  discussed.  First,  the  near 
field  compensation  method  is  explored.  Following  the  discussion  of  this  method,  Matched  Field 
Processing  as  presented  in  [3]  is  discussed,  followed  by  Radial  Reciprocity  as  presented  by 
Kennedy  in  [1,  2],  Robust  Near-Field  Adaptive  Beamfbrming  as  presented  by  Zheng  in  [4],  and 
the  Chirp  Zeta  Transform  (CZT)  as  presented  by  Palmese  in  [5,  6]. 

The  generic  geometry  of  the  test  location  as  well  as  instrumentation  and  signals  used  is 
presented.  More  detail  on  the  test  site  geometry  is  presented  when  each  different  array  and 
source  position  is  discussed  later.  The  array  design  frequency  is  also  presented. 

The  processed  results  of  the  test  data  is  discussed.  Deep  water  results  using  a  far  field 
conventional  beamformer  are  presented  first  to  establish  a  baseline.  The  results  using  the  near 
field  beamforming  algorithm  is  presented  next  to  confirm  that  it  works  property.  The  array  is 
then  moved  into  the  source’s  near  field  and  both  conventional  beamforming  and  near  field 
beamforming  results  are  presented.  The  results  of  the  array  element  localization  calculations  are 
then  explored  and  the  method  of  correction  for  the  error  in  the  vertical  array  assumption  is  given. 
The  results  of  the  deep  water  test  using  the  array  curvature  modification  is  then  presented.  The 
work  then  moves  into  shallow  water  using  the  conventional  beamforming  algorithm,  the  near 
field  algorithm  without  accounting  for  array  curvature,  and  the  near  field  algorithm  with  the 
array  curvature  modification. 

Conclusions  are  drawn.  Future  work  and  possibilities  are  also  discussed. 
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RESULTS 


Beamforming  is  a  method  of  spatial  filtering  which  is  analogous  to  filtering  a  narrowband  signal 
in  the  frequency  domain.  It  is  used  with  an  array  of  sensors  to  accomplish  the  following  goals  [7]: 

1)  Reduce  the  ambient  noise  with  which  a  plane  wave  signal  must  compete. 

2)  Permit  the  separate  detection,  or  resolution,  of  plane  wave  signals  arriving  from  different 
directions. 

3)  Permit  the  measurement  of  the  direction  of  arrival  of  plane  wave  signals. 

Spatial  filtering  is  compared  with  frequency  domain  filtering  in  Figure  2.1. 

Array  bcanifonning  is  used  for  source  detection,  localization,  and  identification.  Many  differ¬ 
ent.  bcaniforming  methods  have  been  developed  over  the  years. 


Far  Field  Beamforiuing 

Far  field  Ixxun  forming  and  line  array  processing  are  explored  in  great  depth  in  many  textl>ooks. 
This  paper  follows  the  method  shown  in  Kinslcr,  Ftey,  Copjxuis,  and  Sanders  [8,  p.  195].  Consider 
a  line  array  of  hydrophones  as  shown  in  Figure  2.2,  with  N  elements  and  with  separation  distance 
d  [8].  If  a  source  generates  a  pressure  wave  of  the  form  ^7  j  exp  [7  (urf  —  JtrJ)]  that  is  received 
by  the  ith  source  (where  r<  is  the  distance  from  the  source  to  the  tth  clement),  the  averaged 
pressure  received  by  the  array  would  be: 


(2.1) 

N  ^  r  g 
1 

Lot  us  define  n  -  r\  —  -  (i  —  1)  Ar  where  Ar  =  d sin  0.  r  represents  an  average  distance 
1  between  the  source  and  the  receiver  elements.  Ar  represents  the  incremental  distance  lx*twocn 
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Figure  2.1:  Comparison  of  frequency  and  spatial  domain  filter  functions:  (a)  power  spectral 
density  of  sine-wave  signals  in  broadband  noise;  (b)  angular  density  of  acoustic  plane  wave 
signals  in  ambient  noise.  (Tbken  from  [7]) 

receiver  elements  toward  the  source.  0  is  the  steering  angle. 

If  the  assumption  is  made  that  the  source  is  in  the  far  field,  where  r  >  (jV  -  l)d  (where 
the  distance  from  the  array  to  the  source  is  much  larger  than  the  length  of  the  array),  then  the 
raypaths  between  the  source  and  each  receiver  are  relatively  parallel.  The  distance  from  the 
source  to  each  element  can  be  defined  as  q  =  n  -  (»  -  l)Ar,  («  =  1, 2, ... ,  N).  It  is  necessary 
to  retain  this  relationship  exactly  in  the  phase,  but  in  the  amplitude  terms,  one  may  use  the  far 
field  approximation  ^  cs  This  gives  a  simplified  equation  taking  these  far  field  assumptions 
into  account: 

p  (r,  0,  t)  =  (2.2) 

iael 

One  of  the  primary  forms  of  bcamforming,  and  the  first  method  presented,  is  called  delay  and 
sum  beamfomiing  and  it  can  he  done  in  both  the  time  and  frequency  domain.  The  time  domain 
method  will  be  presented  first.  This  method  delays  cadi  clement  by  a  specific  time  interval  based 
on  the  difference  in  path  lengths  between  elements,  A r,  which  can  lx?  calculated  from  the  element 


8 


1  2  3  4  &  A  ,V 

- L  -  (A  1U - 


Figure  2.2:  Line  Array  Geometry 


separation  and  the  steering  angle  #o* 


A  r  =  dsin0o 

ns  well  as  the  speed  of  sound  in  the  medium: 


(2.3) 


T  = 


A r 

c 


(2.4) 


By  implementing  a  time  delay  and  summing  over  the  entire  array,  the  steering  angle  that 
receives  the  most  energy  can  be  identified,  revealing  that  to  be  the  angle  of  incidence  and  the 
direction  of  the  source.  This  also  serves  to  help  suppress  background  noise  at  angles  other  than 
the  incident  angle.  This  delay  can  be  added  as  a  phase  delay,  e**1,  where  =  iu)T  =  ifcrfsinflo, 
to  Eq.  2.2  to  get: 


p{r,e,t,9o)  = 


(2.5) 


1=1 


Delay  and  sum  beamforniing  can  also  be  performed  in  the  frequency  domain.  The  equation 
for  this  can  easily  be  found  by  performing  a  Fourier  Transform  on  Eq.  2.5: 
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$  /  N  \ 

F(r,0,u,Oo)  =  ^  I  - -c  JW<e#<wl  *r>r~J(^)*Ar  f  1lfcArH  dt 

~T 

J  dl  (2.G) 


This  is  an  unstcered  Fast  Fourier  Transform  multiplied  by  In  [9],  a  weighting 

factor  is  added  to  the  array  to  control  the  siddol>c5  of  the  beam  pattern. 


AT- 1 

F(M)  =  V  WiXi(k\e-*likAri  (2.7) 

i»0 

where  ir;n  is  a  weighting  coefficient  hikI  A  i(k)  are  the  contplex-vuluecl  coefficients,  or  our  data 
[joints.  Tliis  equation  was  then  vectorized  in  [10],  making  the  implementation  of  this  method 
significantly  easier  and  much  more  computationally  manageable: 


r(/)  =  «lH(0,/)Wx(/)  (2.8) 

Here,  0  is  the  incident  angle,  W  is  a  weighting  matrix: 

W  =  diaR{ir, . \Vn}  (2.9) 

dH(0,  /)  is  a  steering  vector  dependent  on  incident  angle,  frequency,  element  spacing,  and 
spctxl  of  sound: 

dH(<?,/)  =  [<J2*/Ar1<j2ir/Arv]r  (2.10) 

and  x(/)  is  the  data  vector: 

*(/)  =  1*.  (/)...  An(/)f  (2.11) 

Eq.  2.S  can  produce  a  Bearing- Time  Record  (BTR),  which  show's  where  the  signal  exists  in 
space  and  time.  This  can  be  especially  useful  in  detecting  the  movement  of  a  source. 

Figure  2.3  is  a  BTR  which  displays  the  bearing  angle  on  the  x-axis  and  time  on  the  y-axis. 
The  color  axis  here  displays  the  sound  pressure  level  in  dB  re  1  //Pa.  It  can  be  effectively  used  to 
show  vertical  (or  horizontal,  dcjxmding  on  the  setup  of  the  array)  movement  of  sounds  over  time. 
If  the  source  is  not  moving,  a  BTR  can  be  collapsed  into  a  more  easily  readable  two-dimensional 
plot  that  displays  sound  pressure  level  vs.  bearing  angle.  This  is  shown  in  Figure  2.4. 
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Figure  2.3:  Sample  Bearing  Time  Record 
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Figure  2.4:  Sample  Averaged  Bearing  Time  Record 
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Using  these  graphs,  while  the  source  level  can’t  be  known  without  the  range,  the  angle  of 
incidence*  can  be  easily  identified. 


Bartlett  Beamforniing 

Another  very  similar  algorithm  is  presented  in  [3].  This  Bartlett  bcamformcr  takes  the  same 
approach  of  adding  the  shaded  signal  from  each  hydrophone  together  to  find  the  incident  angle 
of  the  signal: 


Bll<xri(Om) 


*0  =  1 


(2.12) 


Here  s,j  +  riij  are  signal  and  noise  that  are  part  of  a  Cross-Spectral  Density  Matrix  and  tv,  is 
the  tth  element  of  a  steering  column  vector.  Eq.  2.12  can  be  vectorized  using  the  steering  column 
vector  w  and  the  Cross-Spectral  Density  Matrix  K.  The  CSDM  can  be  constructed  using  data 
vretnrs  d  and  using  f  to  represent  the  complex  transpose  operation: 


K*  =  (dd+)f  =  dd*  =  K  (2.13) 

Using  these  vectors  and  matrices,  Eq.  2.12  becomes: 

BuartiO')  =  wt(0a)K(<W)w(0.)  =  w*Kw  (2.14) 

Ill  this  aquation,  0irttc  indicates  the  direction  from  which  the  signal  is  incident.  This  equation 
again  provides  a  plot  showing  the  incident  angle  of  the  source. 

Adaptive  Beamforniing 

Adaptive  beaniforining  is  a  method  of  supressing  the  sidtdobes  of  a  traditional  beamfornier.  It 
uses  the  incoming  signal  to  build  the  weight  vectors.  There  are  many  different  types  of  adaptive 
lieamforuiing,  but.  the  one  presented  here  is  the  Minimum  Variance  distortionless  processor  (MV). 

First,  a  weight  vector  needs  to  1m?  defined  that  will  attenuate  the  bcamfonner’s  output  in  every 
direction  except  the  signal’s  incident  angle.  The  weight  vector  that  minimizes  the  functional  can 
l>e  written  as: 


F  =  w|(l,Kw,Wv  +  7(wtwvw  -  1) 


(2.15) 
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The  second  tenn  in  this  equation  utilizes  a  Lagrangian  multiplier  7  in  order  to  indude  the 
unity  gain  in  the  look  direction.  Several  properties  and  methods  can  be  used  to  solve  for  wjifv 
here.  First,  K  is  Hermit ian;  that  is,  it  is  a  square  matrix  that,  is  equal  to  its  own  conjugate 
transpose.  Next,  because  a  Lagrange  multiplier  is  being  used,  the  gradient  can  be  found  with 
respect  to  the  variable  wyr  and  it  can  be  set  equal  to  zero: 


2Kwj \4V  4  7W  =  0  (2. 1C) 

This  gives: 

W MV  =  — ►  W^v  =  --(K-'w)t  (2.17) 

Using  the  second  term  of  Eq.  2.15  as  the  unity  gain  in  the  look  direction: 

=  wjfvw  -1  =  0  (2.18) 

Combining  these  equations  results  in: 


7  =  — 2(w*K-1w)-1 

This  equation  can  then  be  used  to  find  an  expression  for  w \tv: 


WAfV 


K_l  w 
w*K  l\v 


(2.19) 


(2.20) 


Tliis  weight  rector  follows  the  standards  of  adaptive  hcainfonning  in  that  it  utilizes  the 
received  data  to  construct  tlie  weight  vector  that  is  then  applied  to  the  received  data.  This 
wight  vector  can  now  be  used  in  Eq.  244: 


Bmv(6.)  -  [wf(0,)K  ,(0tr«)w(^))  1  (2.21) 

The  Minimum  Variance  distortionless  beaniformcr  is  a  higli-resolntion  beamformcr.  It  is 
designed  to  suppress  the  sidolol)cs  and  have  a  narrower  main  lol>e  than  typical  plane  wave  beam- 
formers.  While  it  performs  this  function  very  well,  it  is  most  effective  with  correlated  noise, 
wliich  is  not  present  in  the  cases  in  this  paper. 


Near  Field  Beainforming 

Far  field  bcamforming,  while  quite  convenient  and  accurate  in  far  field,  free  field  measurements, 
loses  its  precision  in  near  field,  multipath  environments  where  the  far  field  assumptions  cannot  !>e 
made.  Many  different  methods  have  been  developed  to  simplify  the  near  field  problem.  Several 
of  these  will  be  discussed,  starting  with  the  method  chosen  for  tliis  thesis. 
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Near  Field  Compensation 

Tills  method  accounts  for  the  spherical  nature  of  wavefronts  and  calculates  tin?  path  length  dif¬ 
ference  to  each  element  in  the  array.  It  uses  the  same  method  as  the  delay  and  sum  bcainforming 
in  the  frequency  domain  from  the  beginning  of  the  chapter,  but  without  making  the  far  field 
assumptions. 

Starting  with  Eq.  2.1,  far  field  assumptions  cannot  be  made,  so  this  equation  becomes: 

A  1 

p(r,0,t)  =  AciutY,—  ei{*,~ku)  (2.22) 

.-i  r' 

Note  that  a  phase  delay  &  has  been  added  which,  as  in  the  far  field  case,  can  be  defined 
as  =.  *7rAr*.  The  solution  for  A for  the  near  field  case  will  be  presented  later.  Again,  the 
frequency  domain  is  the  domain  of  choice  for  this  thesis,  so  an  FFT  of  tliis  equation  gives 

i  rT/2  /  N  i  \ 

F(u)  -  I  Ac**1  Y  -|rr*)  e~iuidt  (2.23) 

2?r  J  r/2  Vtl  r*  / 

and  c<in  lx*  simplified  to 

1  1  rTP 

F(u9)  =  /  dt  (2.24) 

2x  “  Ti  J-T/2 

Once  again,  it  appears  the  only  addition  to  the  unsteered  FFT  is  a  phase  delay  term  that  is 
now  dependent  on  l>oth  steering  angle  and  range.  To  find  that  range,  the  path  length  difference 
must  l>e  recalculated  for  each  element.  By  varying  lx>th  slant  range  and  look  direction,  the 
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range  to  each  element  can  be  found  using  the  known  array  geometry  (Figure  2,5)  and  the  Law 
of  Cosines 

r{  =  y/r2  +  6?  -  2r6<  cog(90  +  ^,)  (2.25) 

where  6,  =  (i  —  l)d  and  is  the  steering  angle.  Ar*  can  then  lie  computed 

Art  =  r*  —  r  (2.26) 

This  Ar*  can  l>e  inserted  into  Eq.  2.10  and  subsequently  Eq.  2.8.  This  will  compute  a  BTR 
for  a  fixed  range.  By  varying  r  in  Eq.  2.25,  a  map  of  possible  locations  of  the  source  based  on 
slant  range  and  look  direction  is  created  Figure  2.C  shows  the  output  of  this  algorithm. 

The  location  of  the  maximum  of  this  plot  is  the  location  of  the  source  and  is  represented  by 
the  diamond.  The  BTR  can  be  examined  at  the  peak’s  range  if  we  are  interested  in  the  evolution 
of  the  (signal Ts  location  over  time. 

This  is  a  very  simple,  yet  extremely  flexible  method  of  beamforming,  hence  the  reason  it  is 
used  in  this  paper.  Some  other  methods  are  now  presented  to  show  the  different  algorithms  for 
near  field  beamforming  that  have  been  established. 

Matched-Field  Processing 

Matched-Field  Processing  (MFP)  is  a  processing  technique  that  adds  an  additional  dimension  to 
the  traditional  plane  wave  beamformer.  Jensen  covers  the  subject  of  MFP  in  [3].  In  both  cases, 
the  data  is  matched  to  a  replica  field.  In  Conventional  Beamforming,  that  replica  is  plane  waves, 
where  in  Matched-Field  Processing,  the  replica  consists  of  the  theoretical  field  using  dimensions 
of  angle  and  range.  This  allows  the  MFP  to  localize  in  angle  and  range  as  opposed  to  only 
detecting  an  incident  angle,  as  is  done  in  CBF.  Instead  of  only  scanning  through  steering  angles. 
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Figure  2 .G:  Sample  Near  Fidel  Plot 


a  generated  source  is  placed  at  various  points  in  an  angle /range  search  grid.  The  transmission 
path  is  then  onlculaUxl  using  standard  propagation  models,  and  the  tbcorrtical  result  is  placed 
into  a  replica  matrix.  This  matrix  is  then  correlated  with  the  real  data,  and  the  true  location  is 
at  the  peak  of  the  correlation. 

The  MFP  algorithm  can  he  modeled  just  like  a  plane  wave  heamfortner  such  as  in  Eq.  2.14. 
Now,  instead  of  constructing  the  weight  vector,  w,  using  a  simple  far  field  approximation  method, 
the  wwv  equation  can  be  used  to  come  up  with  the  theoretical  results  for  cadi  source  search 
location,  a.  K(a<rt^)  represents  the  rc?al  data  that  is  matched  against  the  replica. 

BuartU1)  =  wt(«)K(»en..)w(a)  (2.27) 

Hero,  the  rnax  of  /}#««-*(**)  the  real  source  location,  a tw  As  in  the  CBF  case,  the 
MFP  possesses  sidelobes  that  are  present  in  the  form  of  ^ambiguous  peaks.”  These  peaks  can 
Ik?  suppressed  by  adaptive  algorithms,  such  as  the  Minimum  VWiaiice  processor,  by  adding  the 
range  dimension  to  the  algorithm. 

tf.uv(n)  =  [w+(n)K~  ‘(a<w)w(a)]  1  (2.28) 

Once  again,  the  real  source  location,  a iruc,  is  located  at  the  max  of 

This  only  scratches  the  surface  of  the  realm  of  MFP,  but  [3]  has  a  list  of  various  applications 
ami  references  for  many  of  the  algorithms  that  have  been  developed. 
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Radial  Reciprocity 

One  method  of  near  field  benmformiiig  that  was  developed  by  [1,  2]  uses  a  near  field- far  field 
radial  beam  pattern  transformation.  This  method  takes  a  near  field  beampattem  and  converts 
it  to  an  equivalent  far  field  beampattem  for  a  specific  frequency*,  but  this  can  be  expanded  to 
broadband  usage. 

They  start  with  a  beampattem  b  —  6(r,  0 ,  <£,  t)  as  a  solution  to  the  wave  equation  expressed 
in  spherical  coordinates: 


1  0  (  jdb\  1  d  f  .  db\  1  {fib  1  (fib 
r2  Or  \  dr)  *  r2sinO  OO  y^Od)  +  r2sin20  0<^  ~  c2  Ot2 

The  solution  is  then  presented  as: 

br(0, 4>:  k)  =  r-VM  ( ^  anH™l/J.kr)Pn( cos0)+  . . . 

Vn^-0 

OO  OC 

E  E  x  +  W**) 

n— 1  m  =  1 


(2.29) 


(2.30) 


when'  k  =  —iC  js  the  wavenumber,  Pn(*)  is  the  nth  order  Legendre  function,  P™(*)  are  the 
associated  Lrgendrr  functions,  and  the  half  integer  order  Hankel  function  of  the  first  kind  is 
written  as: 

rfl! i./i(*)-^+./iO+ir«+i/3(-) 

where  Jn+1^(*)  is  a  half  integer  order  Bessel  function  of  the  first  kind,  and  j(-)  is  a  half 

integer  order  Neumann  function  (also  a  Bessel  function  of  the  second  kind). 

The  coefficients  an,  /?nm,  and  7nm  can  be  defined  as: 


n-f « 


—  X  f  f  br(6,$;k)Pn  cos  0)  tin  6  dG  d$ 

J  o  J o 


= 


(2.31) 


2xr(«/2)//0)i/2(A-r) 

0nm  = - ; - -1 - ^  x  T  br(0,  k)P?{c<x  0)  sin  0c~int*  dO  d<*>  (2.32) 

2xr~^nm^lxn(kr)  (n  +  m)!  J0  J0 

- +J - -_n — m)  f7T  f  b  ^ ^.k)P™(c<Xi0)ain0cim*d9d4>  (2.33) 

2jrr-0/2)f/(‘)1/2(Jtr)(n  +  n.)!  J0  J0 


Tnrn 


Because  the  coefficients  in  Eq.  2.30  fully  dcseriln?  the  beam  pattern  at  any  range,  it  can  be 
rewritten  for  sources  at  any  point  in  space,  allowing  a  near  field  l>campattcm  to  be  transformed 
to  a  far  field  beampattem,  moving  from  a  finite  sphere  in  the  near  field  to  an  infinite  sphere  in 
the  fcur  field. 

The  following  steps  outline  tlie  procedure  of  using  this  transformation. 

1)  Find  ow,  /?wm,  and  7mn  for  the  wanted  nearficld  beampattem  brd(0^;k)  where  r  =  rrf. 

2)  Use  Eq.  2.30  to  calculate  b^d^;  k)  at  r  =  oo. 
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3)  Di’sign  a  fnrfield  iK'ninforincr  to  create  this  desired  lieauipattern. 
The  far  field  synthesis  equation  is  established  in  [1]: 


boo(0,<t>;k)  ~  >  (  £  an(-j)n+lrn(cosO)  +  . .  (2.34) 

A  special  case  for  linear  arrays  is  also  explored.  This  rrduces  Eqs.  2.30  and  2.31  coefficients 
to  two  very  simple  equations: 


br(0;k)  =  £>„>•  '/2H{n'ly2(kr)Pn(ctx,6) 


(2.35) 


*=0 


n  -f 


—  f  bt.(0;k)Pn(co$6)sm0d9  (2.3G) 

fcr)  Jo 


This  results  in  a  much  simpler  synthesis  equation: 


bac(6;k)  ~  t  fc°>  J2  “n(-j)n+lrn(<™6)  (2.37) 

n*0 

This  gives  a  simplified  method  of  creating  a  radial  bcampattcni  near  Bold-far  field  transfor¬ 
mation  for  a  specific  case.  Further  details  can  l>e  found  in  [1,  2]. 

Robust  Noar-Field  Adaptive  Beamforming 

Near  field  adaptive  liearnforming  was  explored  by  Zheng,  Gonbran,  and  El-Tanany  in  [4],  They 
tackle  the  robustness  problem  of  creating  an  adaptive  lx^amforaicr  for  the  near  field.  They  go 
on  to  cite  the  advantage  of  distance  discrimination  using  their  algorithm,  saying  uThe  proposed 
l>eam former  can  simultaneously  improve  the  distance  discrimination  and  the  robustness  against 
location  errors  without  additional  constraints.” 

Zheng  starts  by  defining  a  boamformer  with  elements  located  at  xm  =  (rw,0m,<£m)  in  the 
spherical  coordinate  system,  as  shown  in  Figure  2.7,  where  m  =  1,2, ...  t  A /  and  where  M  is  the 
number  of  elements  in  the  array.  Note  that  the  meanings  of  0  and  <p  represent  the  azimuthal 
angle  and  polar  angle  respectively.  Those  definitions  are  common  in  mathematics,  however  they 
arc  opposite  the  typical  convention  found  in  physics.  This  lias  no  bearing  on  the  equations 
ill  this  section.  They  also  define  K  as  the  number  of  taps.  Assume  the  source  is  located  at 
x*  =  (i%,  and  that  it  is  within  onr  near  field  parameters  as  defined  in  Eq„  1.1.  The  near 

field  steering  vector  would  need  to  lie  used 


r/2x/r./r 


/ri,/c  p/2»/(rm*/c—  fr)  j./c-K-fl) 


n- 


rxi* 


where  T  indicates  the  transpose  operation,  /  is  the  operating  frequency',  c  is  the  sound  speed. 
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Figure  2.7:  Constrained  spatial  region  fi  =  (±Ar,  ±A0,  d:A<£).  The  presumed  focal  point  is 
xp  =  (rpyOp^p).  (Figure  taken  from  [4]) 

r„  =  |x*|  is  the  absolute  distance  from  the  source  to  the  center  of  the  array,  and  r™  =  |xm  —  x*| 
is  the  absolute  distance  from  the  source  to  the  mth  dement  of  the  array.  If  the  input  to  the 
beamformer  is  allowed  to  be  samples  in  a  vector  N  =  MK  and  define  it  as  u(fc),  where  k  is  the 
tap  index.  The  output  of  the  beamformer  can  be  defined  simply  as  y(fc)  =  w^u(fc),  where  w  is 
the  weight  vector  and  H  is  the  complex  conjugate  transpose  operation.  The  linearly  constrained 
minimum  variance  (LGMV)  beamforming  method  can  then  lie  used  to  minimize  the  output  power 
around  constraints 


min  [w"Rouw]  (2.39) 

subject  to  CWw  =  li  (2.40) 

where  R^u  is  the  N xN  covariance  matrix  of  the  input  vector,  u.  In  this  case,  h  is  the  wanted 
response  vector  and  C  is  the  matrix  of  constraints.  This  gives  an  optimal  solution 

=  R-»-,C(C<rRM(1-IC)-1h  (2.41) 

which  can  be  broken  into  two  components 
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w,  =  C(C,/C)"1h 

wa  =  (C"RuuCn)-‘C"Rullw7 


(2.42) 

(2.43) 


where  \vq  is  the  fixed  bcamfonner  component,  \va  is  the  unconstrained  adaptive  component, 
and  Ca  is  the  signal  blocking  matrix.  The  solution  to  wft  using  the  normalized  least  mean  square 
(XLMS)  algorithm  can  be  written 

w-(* + 1}  =  w*(fc) + (244) 

where  the  error  is  e(fc)  =  [w^  —  Cawa]*  ti(fc)  and  the  step  size  is  fi.  The  robustness  issue 
can  now  l>e  exploml.  Consider  a  signal  in  a  region  fi  =  (rfcAr,  :tA#,  ±A$)  arenmd  a  focal  point 
xy  —  (rj.%  0p,  and  a  frequency  hand  B .  If  the  source  sample  vector  is  defined  as  s(fc)  and 
the  power  spectrum  density  as  S(/),  then  the  covariance  matrix  can  be  written  as 

R«*  =  E{s{k)sT(k)}  =  11  I  I  5(/)a(x./)a"(x,/)dfrfx  (2.45) 

where  a(x,/)  is  defined  in  Eq.  2.38.  A  Karhnnen-Loew  expansion  vatu  be  performed  on  the 
sample  vector  s(k)  and  becomes  s(fc)  —  ^2r^~\  snvn  where  sn  =  v^s(Ar).  Approximating  this 
expansion  gives 


L 

s(k)  m  s (k)  =  s„v„  (2.4G) 

n=l 

The  eigenvalues  An  correspond  to  the  energy  of  the  sample  vector  projected  onto  vn.  L  is 
chosen  to  obtain  an  approximation  error  that  can  be  written  as 

<{L)  =  E{\s(k)-s(k) |2}  =  Y,  vpt»v„=  Y  A"  (247) 

The  eigenvalues  of  a  broadband  signal  doorcase  quickly  when  n  <  fi,  where  0  =  2R7’(xB), 
and  T(xJ  is  tlie  length  of  the  time  window  as  a  function  of  the  source  location 

r(x8)  =  r(x«)  +  (2.48) 

where  r(xB)  =  [max(rmjw)  —  min( rm9)\fc/FB  where  Fm  is  the  sampling  frequent  The  con- 
strain!  design  method  can  now  be  established.  If  a  flat  spectrum  in  B  is  assumed,  can  be 
easily  solved  for 

\  1  J  1 

R~  =pEE  «(^/j)«"(^./j)  =  pAAr  (2.49) 

f-l 

wliere  P  =  /  J  and  A  is  defined  as 


20 


j\  [nc(xj ,  f\ ),***,  <V (Xj ,  fj * v{^/ ,  fj )  j /l  )i  •  •  •  i  /jr ^  1  //  )j  { 2 . ‘->0 ) 

where  a c(x,,  /,)  and  a,(x#1  /j)  are  tlie  real  and  imaginary  parts  of  the  steering  vector.  Mean¬ 
while  the  constraint  equation  can  be  written  as 


Arw  =  g 


(2.51) 


where  g  Ls  the  response  vector 


g  =  \gn  ros(27r/,7-,), . . . ,  g,j  co^{2~fjT,)\gu  sin(2^/,r,), . . .  ,fl/jsin(2Jr/Jr,)]r  (2.52) 

r,  is  the  group  delay  of  the  signal  at  location  x*  relative  to  the  origin  and  gij  arc  the  desired 
amplitude  responses.  The  singular  value  decomposition  (SVD)  of  A  can  be  used  to  find  the  low 
rank  representation  of  It**, 


A  =  V£Ur  (2.53) 

where  £  is  a  diagonal  matrix  of  the  values  cr],<X2, . . .  of  A  in  a  descending  order  and  the 
columns  of  V  and  U  are  the  singular  vectors.  The  values  of  A  are  related  by  A„  =(T'//>  to  the 
eigenvalues  of  Therefore, 


A  «  A L  =  VLVLUl  (2.54) 

where  the  columns  of  and  U/,  are  the  L  columns  of  V  and  U  for  the  L  singular  values 
and  —  diag{ai.a2,...,<rt}-  Combining  Eqs.  2.51  and  2.51  gives  V£w  =?  which 

allows  a  constraint  expiation  to  l>e  designed 

C  =  VL 

h  =  ^,Ul'R  (2.55) 

Tlie  procedures  to  design  this  beamformer  are  summarized  by  |4]  as  follows: 

1 )  Select  a  large  number  of  points  /  mid  J  in  a  specified  constraint  region  Q  and  frequency 
hand  B ,  respectively;  Form  the  matrix  A  as  in  2.50  and  the  desired  response  vector  g  as  in  2.52; 

2)  Perform  the  singular  value  decomposition  (SVD)  of  A.  Find  the  matrices  V,  £,  and  U  as 
in  2.53. 

3)  Specify  {in  approximation  error  c(L).  For  example,  an  error  level  of  -40  dB  is  sufficient 
for  unit  gain  constraints.  Then  find  the  required  number  of  constraints  L  by  computing  the 
approximation  error  as  in  2.47; 
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4)  Form  St ,  Vt>  and  Ut  by  choosing  the  largest  L  singular  values  in  S  and  the  corresponding 
L  column  vectors  in  V  and  U,  respectively; 

5)  Compute  the  linear  constraints  C  and  li  as  in  2.55; 

G)  To  implement  the  adaptive  ljeamfonner  in  the  GSC  structure,  compute  \vq  as  in  2.43  and 
wa  as  in  2.43  or  2.44. 

This  method  of  near  field  l»eajnforming  is  based  on  the  far  field  LCMV  beainforcning  al¬ 
gorithm.  While  tliis  seems  to  l>e  an  effective  method  of  near  field  bcamforming,  the  LGMV 
l>eainfomicr  with  Zheng’s  robustness  algorithm  is  not  the  method  chosen  for  this  thesis. 

Chirp  Zcta  Transform 

Another  method  of  near  field  3-D  imaging  exploits  the  computational  ease  of  the  Chirp  Zeta 
transform  (CZT)  algorithm  and  was  presented  in  [5]  and  [6],  A  two  dimensional  army  must  be 
used,  as  scanning  in  tlirco  dimensional  space  requires  two  different  steering  angles.  Thus,  a  planar 
array  composed  of  M  x  jV  elements  is  used.  The  previously  defined  d  is  the  separation  distance 
l)etweeii  elements.  Following  time-domain  l>eamfonning  convention  for  two  steering  directions, 
the  beam  signal  can  l>e  written  as 


Af  —  1  S~\ 

b(t,0ap,0cq)  =  ^  Wm,nSm.w(t  ~  ^(cq,  0ap ,  6eq,  m,  Tl))  (2  LG) 

m=0  n^O 

for  0  <  p  <  Mb  —  1  and  0  <  q  <  Nt  -  1  where  p  and  q  are  the  beam  signal  indexes  and 
sulweripts  a  :uid  c  refer  to  the  azimuth  and  elevation  directions.  Also,  the  time  delay,  r,  can  lie 
expressed  as 

T(ro,0ap,Otq,m,n)  —  (2.57) 

where 


X  =  r5  +  ((m  -  M/2)df  4-  ((n  -  N/2)d)2  -  2 (m  -  M/2)ro<f  sin  0ap  -  2(n  -  N/2)r0dsin  0C, 

and  sm.n( 0  is  Ule  incoming  signal,  wltile  irm>n  is  tbe  weight  factor  that  controls  the  sidelolx? 
levels.  As  usual,  c  is  the  sound  speed.  If  the  Fresuel  approximation  is  assumed,  r  becomes 


T(r0, 


(m  —  M/2)dsiuOap 


-f 


(n  -  Ar/2)d sin  6cq  ((m  -  M/2)df  -  ((n  -  Ar/2)d)2 

c  2rpc 


(2.58) 


The  signal  Is  then  broken  into  snapshots  of  length  I\  and  converted  into  the  frequency  domain. 
This  results  in  the  l>eam  signal  lf{t-,Oap,0C9)  becoming 
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V—.  /  (mdsmOap  ndsmOc„\\ 

D{k,eap,Ocq)  =  y'  Y]  ivtn,nSln,n(k)  •  exp  I  -]2nfk  ( - -  + - —  )  )  • 

n,-0n»0  V  \  C  C  JJ 

■  cxp(j2»  /i-7(m,  n,p.q))  (2.59) 


where 


7K n  p  q)  =  (jmd* mO.,y  4  (JV/2Mg<  +  ((*■  Ji//2)d)2  +  ((n  -  V/2)rf)3  2 ^ 


2roc 


croc 


Not c  that  in  Eq.  2.G0,  the  first  two  terms  are  dependent  on  steering  angles,  hut  not  on  the 
indexes.  For  now  those  terms  can  he  thrown  out  and  revisited  later.  This  now  gives  a  simplified 
fonn  of  Eq.  2.G0 


,  ,  ((m-,\//2)«f)2  ,  ((»*  —  AT/2)rf)2 

7  (m,n)  = - - - + 


2ror 

Eq.  2.59  can  lx*  further  simplified  by  defining 


2ror 


(2.C1) 


(  n  /  dsin  0Up 

=  «P  I -j2*A - 

t Yr*.n  ~  CXp(jf23rA^Vn,n  ) 


(2.C2) 

(2.C3) 


and  rewriting  it  as 


,U-1  AT-1 


*  (*. 0a„. 0rq)  =  53  53  (2.&4) 

nj  =0  n  t) 

If  fc  is  held  constant,  the  CZT  algorithm  can  1x3  utilized.  This  allows  a  scries  of  variables  to 
lx?  set  relating  to  the  azimuth  and  elevation  angles: 


A*  —  exp  j2 xfk-  sin^di^ 

iV„  -  exp  ^j2x/fc^As„  j 


(2.C5) 


As&  — 


si  ndaf  —  sin  0ai 


and 


Ae  =  exp  (~j2xfk~  smO„^ 
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(2.GG) 


He  =  exp  ^j'2« 

A  _  andtf -sinOei 

a*c  “  \rb  i 

where  the  i  and  /  additions  to  the  subscripts  indicate  the  initial  and  final  steering  angles. 
These  angles  can  bo  used  to  find  the  angles  6ap  and  Ocq 


sin  6ap  =  sin0„i  -f-  A sap 
sin0cf  —  sinfle,  4*  A srq 


p  =  0, . . . ,  Mb  —  1 
q  =  0, . . . ,  Nb  -  1 


Finally,  Eq.  2.59  can  be  written  in  a  form  using  the  CZT  algorithm 


(2.07) 

(2.G8) 


a  *Vf — 1 AT—  1  _mj  ^3  ^  »*»)2  ^ 

B{k,oap.orq) «  £  I]  *w  wf5" wr^~ irr^  (2.09) 


*W-1  iV-1 

EL 

m=4)  n  *t) 

Tliis  equation  can  be  represented  by  a  convolution  of  two  matrices 


C(k)  =  vm.KSm.n(k)A?A:W^-Wr£  (2.70) 

D  = 

This  allows  for  fast  computing  using  convolution  methods.  After  nil  the  liearas  are  calculated, 
an  inverse  Fourier  Transform  will  reveal  the  l>enm  signal. 

This  method  is  effective  and  is  a  computationally  easy  load.  However,  this  study  involves  data 
collected  using  a  1-D  array,  and  this  algorithm  is  designed  for  a  planar  array.  Tliis,  therefore,  is 
not  the  method  that  was  chosen  for  this  work. 

Many  different  post-processing  methods  were  presented  in  this  chapter.  Of  the  methods 
presented,  the  ones  chosen  for  this  work  were  the  basic  frequency  domain,  far  field  benmfonning 
and  the  near  field  compensation  method  of  near  field  bcamfonning.  These  will  be  used  in  Chapter 
4  to  analyze  the  data  that  was  collected.  The  next  chapter  will  describe  the  instrumentation 
use nl  and  the  gminctrirs  of  the  test  layouts. 
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IMPACT 

This  study  examined  data  acquired  at  the  Navy’s  Acoustic  Research  Detachment  at  Lake  Pend 
Oreille  in  Bayview,  Idaho  and  processed  the  data  using  a  near  field  compensation  beamforming 
algorithm  with  a  modification  to  account  for  array  curvature.  The  improvements  in  the  accuracy 
of  this  modification  were  presented  as  well. 

Several  methods  of  both  far  field  and  near  field  beamforming  were  presented,  including  the 
methods  used.  Each  method  has  its  benefits  as  well  as  its  disadvantages,  and  the  simplest 
method  with  the  ability  to  adjust  to  a  non-linear  array  was  chosen. 

The  array  and  test  site  geometry,  as  well  as  the  various  equipment  that  was  used  in  testing,  was 
briefly  discussed.  The  many  signals  that  were  transmitted  were  also  presented. 

Results  of  the  processed  data  were  presented.  Results  were  first  presented  for  the  case  of  an 
array  in  deep  water  with  the  source  in  the  far  field.  A  method  of  far  field  beamforming  was  used 
to  determine  the  location  of  the  source,  followed  by  results  of  using  the  near  field  compensation 
method  of  beamforming  to  be  sure  the  near  field  algorithm  worked  properly  while  the  source 
was  in  the  far  field.  The  array  was  then  moved  into  the  near  field  of  the  source  and  it  was  proven 
that  a  far  field  beamforming  method  would  not  work.  The  near  field  algorithm  was  then  used, 
accompanied  by  a  method  for  accounting  for  a  non-linear  array,  followed  by  a  discussion  of  the 
accuracy  of  this  modification.  The  array  and  source  were  then  moved  into  shallow  water  where 
many  reflections  were  present.  The  near  field  algorithm  with  the  array  curvature  modification 
was  used  to  analyze  a  CW  pulse,  a  LFM  pulse,  and  broadband  noise.  The  results  presented 
showed  that  this  method  improved  the  accuracy  of  the  near  field  compensation  beamformer  for 
all  cases  of  a  CW  pulse  and  LFM  pulse.  Broadband  noise  was  less  accurate,  however  it  was 
postulated  that  this  was  due  to  movement  of  the  array  during  the  course  of  the  data  collection 
event. 
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